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Motivated by subdiffusive motion of bio-molecules observed in living cells we study the stochastic 
properties of a non-Brownian particle whose motion is governed by either fractional Brownian motion 
or the fractional Langevin equation and restricted to a finite domain. We investigate by analytic 
calculations and simulations how time-averaged observables (e.g., the time averaged mean squared 
displacement and displacement correlation) are affected by spatial confinement and dimensionality. 
In particular we study the degree of weak ergodicity breaking and scatter between different single 
trajectories for this confined motion in the subdiffusive domain. The general trend is that deviations 
from ergodicity are decreased with decreasing size of the movement volume, and with increasing 
dimensionality. We define the displacement correlation function and find that this quantity shows 
distinct features for fractional Brownian motion, fractional Langevin equation, and continuous time 
subdiffusion, such that it appears an efficient measure to distinguish these different processes based 
on single particle trajectory data. 



I. INTRODUCTION 

Anomalous diffusion denominates deviations from the 
regular linear growth of the mean squared displacement 
(r 2 (i)} ~ Kit as a function of time t, where the propor- 
tionality factor K\ is the diffusion constant of dimension 
cm 2 /sec. Often these deviations are of power-law form, 
and in this case the mean squared displacement in d di- 
mensions 



(r 2 (t)> = 



2dK a 
T{l + a) 



(1.1) 



describes subdiffusion when the anomalous diffusion ex- 
ponent is in the range < a < 1, and superdiffusion for 
a > 1 [J- The generalized diffusion constant K a has di- 
mension cm 2 /sec Q . We here concentrate on subdiffusion 
phenomena. Power-law mean squared displacements of 
the form (|1.1|) with < a < 1 have been observed in a 
multitude of systems such as amorphous semiconductors 
0, subsurface tracer dispersion Q, or in financial mar- 
ket dynamics [3] ■ Subdiffusion is quite abundant in small 
systems. These include the motion of small probe beads 
in actin networks local dynamics in polymer melts 
@, or the motion of particles in colloidal glasses Q- In 
vivo, crowding-induced subdiffusion has been reported 
for RNA motion in E. coli cells Q , the diffusion of lipid 
granules embedded in the c ytop lasm [13] , the propa- 
gation of virus shells in cells [Tj| , the motion of telomeres 



in mammalian cells 12 1 , as well as the diffusion of mem- 



brane proteins and of dextranc probes in HcLa cells 13]. 

While the mean squared displacement Eq. (|1.1[) is com- 
monly used to classify a process as subdiffusive, it does 
not provide any information on the physical mechanism 
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underlying this subdiffusion. In fact there are several 
pathways along which this subdiffusion may emerge. The 
most common are: 

(i) The continuous time random walk (CTRW) 0] in 
which the step length has a finite variance (<5r 2 ) of jump 
lengths, but the waiting time r elapsing between succes- 
sive jumps is distributed as a power law ip(r) ~ t§ /t 1+ol , 
with < a < 1. The diverging characteristic waiting 
time gives rise to subdiffusion of the form (|1.1[) [3 . The 
subdiffusive CTRW in the diffusion limit is equivalent 
to the fractional Fokker-Planck equation, that directly 
shows the long-ranged memory intrinsic to the process 
[l[. Waiting times of the form ip(r) were, for instance, 
observed in the motion of tracer beads in an actin net- 
work Q. We note that recently subdiffusion was also 
demonstrated in a coupled CTRW model [13] . 

(ii) A random walk on a fractal support meets bot- 
tlenecks and dead ends on all scales and is subdiffusive. 
The resulting subdiffusion is also of the form (|1.1[) . and 
the anomalous diffusion exponent is related to the fractal 
and spectral dimensions, df and d s , characteristic of the 
fractal, through a = d s /df [To| . A typical example is the 
subdiffusion on a percolation cluster near criticality that 
was actually verified experimentally [ToT ]. 

(iii) Fractional Brownian motion (FBM) and the frac- 
tional Langevin equation (FLE) that will be in the fo- 
cus of this work and will be defined in Sec. [TTJ The un- 
derstanding of these types of stochastic motions is up 
to date somewhat fragmentary. Thus the first passage 
behavior of FBM is known analytically only in one di- 
mension on a semi-infinite domain [l7j ; the escape from 
potential wells in the framework of FLE has been stud- 
ied analytically [ill, [l9j] and numerically (2(| [2l| ; and a 
priori unexpected critical exponents have been identified 
for the FLE [22j]. Here we address a fundamental ques- 
tion related to FBM and FLE. Namely, what is their 
behavior under confinement? Two main aspects of this 
question will be addressed. One is the study of the re- 
laxation towards stationarity in a finite box by means of 
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the ensemble averaged mean squared displacement. We 
also investigate a new quantity used to characterize the 
motion, the displacement correlation function. For these 
aspects we also study the dependence on the dimension- 
ality of the motion. 

The second aspect concerns how FBM and FLE mo- 
tions under confinement behave with respect to ergodic- 
ity. Experimentally the recording of single particle tra- 
jectories has become a standard tool, producing time se- 
ries of data that are then analyzed by time rather than 
ensemble averages. For subdiffusion processes both are 
not necessarily identical. In fact for CTRW subdiffusion 
with an ensemble averaged mean squared displacement 
of the form the time averaged mean squared dis- 

placement 

W^T) = f A [x(t + A) - x(t)] 2 dt (1.2) 

on average scales like (s 2 (A,T)^ ~ A/T 1 ' 01 0, H]. 

That is, the anomalous diffusion is manifested only in 
the dependence on the overall measurement time T and 
not in the lag time A, that defines a window swept along 
the time series. Thus time and ensemble averages are 
indeed different. In contrast, for normal Brownian dif- 
fusion {^(5 2 (A,T)^ ~ A is independent of T, and time 

and ensemble averages become identical, i.e., the sys- 
tem is ergodic. Different from CTRW subdiffusion, sys- 
tems governed by FBM or FLE arc ergodic. The er- 
godicity breaking parameter measured from time aver- 
aged mean squared displacements converges algebraically 
to zero [ergodic behavior] as the measurement time in- 
creases, the convergence speed depending on the Hurst 
exponent H = a/2 [25[. For the FLE case, however, it 
was also shown that the crgodicity measured from the ve- 
locity variance can be broken for a class of colored noises 

One of the open questions in the context of ergodic- 
ity breaking for FBM and FLE in the above sense is the 
influence of boundary conditions on the time averages. 
It was shown for CTRW subdiffusion that confinement 
changes the short time scaling ^<5 2 (A,T)^ ~ A/T 1 ~ a to 

the long time behavior (s 2 (A,T)^ ~ (A/T) 1_Q ^Hl- 

Although we expect that FBM and FLE processes be- 
come stationary under confinement and, for instance, 
attain the same long-time mean squared displacement 
dictated by the size of the confinement volume, we in- 
vestigate how fast this relaxation actually is, and how 
it depends on the volume and the dimensionality. To 
this end we study the ergodicity breaking parameter for 
the system. We find that for both FBM and FLE con- 
finement actually decreases the value of the ergodicity 
breaking parameter with respect to unbounded motion, 
i.e., the process becomes more ergodic. Ergodicity is also 
enhanced with increasing dimensionality. We also discuss 
how FBM and FLE motions can be distinguished from 
time series from single particle trajectories. 



The paper is organized as follows. In Sec. UU we in- 
troduce FBM and FLE motions, and review briefly their 
basic statistical properties. In Sec. IIIIl we describe the 
numerical scheme for simulating FBM and FLE in con- 
fined space. Simulations results are presented in Sec. IIVI 
and |Vl where we discuss the effects of confinement and 
dimensionality on time-averaged mean squared displace- 
ment trajectory, ergodicity, and displacement correlation. 
We draw our Conclusions in Sec. IVII 

II. THEORETICAL MODEL 

We here define FBM and FLE. These two stochas- 
tic models share many common features, however, their 
physical nature is different. In the following we will see, 
in particular, how the two can be distinguished on the 
basis of experimental or simulations data. 

A. Fractional Brownian motion 

FBM was originally introduced by Kolmogorov in 1940 
[29^ and further studied by Yaglom |3(|. In a different 
context it was introduced by Mandelbrot in 1965 [3l| 
and fully described by Mandelbrot and van Ness in 1968 
in terms of a stochastic integral representation [32j . In 
the latter reference the authors wrote that "We believe 
FBMs do provide useful models for a host of natural time 
series". This study was motivated by Hurst's analysis of 
annual river discharges (33j , the observation that in eco- 
nomic time series cycles of all orders of magnitude oc- 
cur [34j], and that many experimental studies exhibit the 
now famed 1/ / noise. FBM by now is widely used across 
fields. Among many others FBM has been identified as 
the underlying stochastic process of the subdiffusion of 
large molecules in biological cells [H, Hi, • We note 
that FBM is neither a semimartingale nor a Markov pro- 
cess, which makes it quite intricate to study with the 
tools of stochastic calculus [37| H|[ • 

FBM, x H (t), is a Gaussian process with stationary in- 
crements which satisfies the following statistical proper- 
ties: the process is symmetric, 

(x H (t)) = 0, (2.1) 

with x H (0) = 0; and the second moment scales like 
Eq. fTI) : 

(x H {t) 2 } = 2K H i 2H . (2.2) 

For easier comparison with other literature we introduced 
the Hurst exponent H that is related to the anomalous 
diffusion exponent via H = a/2. The Hurst exponent 
may vary in the range < H < 1, such that FBM de- 
scribes both subdiffusion (0 < H < 1/2) and supcrdiffu- 
sion (1/2 < if < 1). The limits H = 1/2 and H = 1 cor- 
respond to Brownian and ballistic motion, respectively. 
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Finally, the two point correlation behaves as 

{x H (h)x H (t 2 )} = K H {t\ H + t\ H - |ti - t 2 \ 2H ). (2.3) 

Here (•) represents the ensemble average. It is con- 
venient to introduce fractional Gaussian noise (FGN), 
£ H (*), from which the FBM is generated by 



x H (t) 



dt'£ H (t'). 



FGN has the properties of zero mean 

(Z H (t)) = o 

and autocorrelation [39|, 

(Z H (h)$ H (t 2 )) = 2K H H(2H - I))*! - t 2 \ 



(2.4) 



(2.5) 



2H-2 



AK H H\h - t. 



\2H-1 



S(h-t 2 ), 



(2.6) 



as can be seen by differentiation of Eq. (|2.3p with re- 
spect to t\ and *2- Here we see that Kh plays the role 
of a noise strength. For subdiffusion (0 < H < 1/2) the 
autocorrelation is negative for t\ ^ *?, i.e., the process 
is anti-correlated or antipersistent [411 ] . In contrast, for 
1/2 < H < 1, the noise is positively correlated (persis- 
tent) and the motion becomes supcrdiffusivc. For nor- 
mal diffusion (H = 1/2), the noise is uncorrelated, i.e., 
(Z H (h)€ H fa)) = 2K H 5{t 1 -t 2 ). For further details com- 
pare the discussions in Ref. [12, 5H . 

We define c?-dimensional FBM as a superposition of 
independent FBMs for each Cartesian coordinate, such 
that 



**(t) 



d 

E 

i=l 



dfff (f)*i, 



(2.7) 



where Xj is the Cartesian coordinate of the ith component 
and £^ is FGN which satisfies 



<£f (*)> = o 



(2.8) 



and 



<£?(*i)£f(*2)> = 2K H H(2H - 1)1*! - * 2 



2H-2 



5 « 



+ 4A- H if|ti-t2H-M(t 1 -* 2 )«y. 

(2.9) 

From this definition, d-dimensional FBM x H (f) has the 
properties of zero mean 



(x ff (*))=0, 



variance 



(2.10) 



(2.11) 



(x"(*f)=2^ ff *^, 
and autocorrelation 
(x H (t 1 )-x H (t 2 ))=dK H (t 2H + t 2H -\t 1 -t 2 \ 2H ). (2.12) 



Note that \x. H (t)| cannot satisfy these properties, thus it 
is not an FBM. 

A few remarks on this multidimensional extension of 
FBM are in order. We note that, albeit intuitive due 
to the Gaussian nature of FBM, this multidimensional 
extension is not necessarily unique. In mathematical lit- 
erature higher dimensional FBM in the above sense was 
defined in Rcfs. [43, [44[. In physics literature an analo- 
gous extension to higher dimensions was used in Ref. [l3| 
based on the Weierstrass-Mandelbrot method. To ver- 
ify that this d-dimensional extension is meaningful, we 
checked from our simulations of d-dimensional FBM that 
the fractal dimension of FBM, df = 1/H for H > 1/d 
[H| , is preserved in higher dimensions. Moreover we used 
an alternative method to create FBM in (^-dimensions, 
namely, to use ID FBM to choose the length of a radius 
and then choose the space angle randomly. The results 
were equivalent to the above definition to use indepen- 
dent FBMs for every Cartesian coordinate. We are there- 
fore confident that our definition of FBM in d dimensions 
is a proper extension of regular ID FBM. 



B. Fractional Langevin equation motion 

An alternative approach to Brownian motion is based 
on the Langevin equation j46l - |48| 



d 2 y(t) _ _dy(t) 
1 dt 2 



-7" 



dt 



(2.13) 



where £(*) corresponds to white Gaussian noise. 

When the random noise £(*) is non-white, the resulting 
motion is described by the generalized Langevin equation 
(GLE) 



d 2 y(t) 
1 dt 2 



Jfr(t-t')^dt' + m, (2-14) 



where m is the test particle mass, is the memory 
kernel (49l - [5l| which satisfies the fluctuation-dissipation 
theorem (£(*)£(*')) = k B T^Jf(t-t'). When g is the FGN 
introduced above, Jfr decays algebraically, and Eq. (|2.14[) 
becomes the fractional Langevin equation (FLE) 

= -jT(2H - V-^gV® + rt H (t). (2.15) 

Here, 7 is a generalized friction coefficient. We also define 
the coupling constant 



2K H H(2H - 1) 



(2.16) 



imposed by the fluctuation dissipation theorem, and the 
Caputo fractional derivative (52| 



-]2-2H 



-y(t) 



dt 2-2H^ ' T(2H-1)J 



r dt'(t-t') 2H - 2 ^- (2.17) 
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Note that in Eq. (|2.15[) the memory integral diverges for 
H smaller than 1/2, such that the Hurst exponent in the 
FLE is restricted to the range 1/2 < H < 1. 

It can be shown that the relaxation dynamics governed 
by the FLE ([27X5]) follows the form 

(y(t)} = v tE 27l 2 (- 7 t 2l? ) (2-18) 

for the first moment, where vq is the initial particle veloc- 
ity. The rescaled friction coefficient is 7 = 7T(2iJ — 1) / m. 
The coordinate variance behaves as 

{y i(t))=2^E 27l3 (- 1 t^) (2.19) 

where (v 2 ,) = fcgT/m is assumed, and we employed the 
generalized Mittag-Lcffier function [53j 

00 „ 

Ea,p(z) = £ TV \ R , ( 2 - 2 °) 

whose asymptotic behavior for large z is 

n—l v ' 

Thus the mean squared displacement shows a turnover 
from short time ballistic motion to long time anomalous 
diffusion of the form (HD, [55| 

<fW)~{%i,\Zl- (2 ' 22) 

Therefore, for persistent noise with 1/2 < H < 1 the re- 
sulting motion is in fact sub-diffusive, i.e., the persistence 
of the noise has the opposite effect than in FBM. 

In analogy to our discussion of FBM in a d-dimcnsional 
embedding the FLE is generalized to 

where y(t) = £;2/*(*)£i, = an d X is 

the memory tensor which is in diagonal form (i.e., = 
J€(t— t')Sij) in the absence of motional coupling between 
different coordinates. 



III. SIMULATIONS SCHEME 

We here briefly review the simulations scheme used to 
produce time series for FBM and FLE motion. 

A. Fractional Brownian motion 

d-dimensional FBM is simulated via Eq. ()2.7|) by nu- 
merical integration of ^(i). The underlying FGN was 



generated by the Hosking method which is known to be 
an exact but time-consuming algorithm [56^ . We checked 
that in the one-dimensional case the generated FBM in 
free space successfully reproduces the theoretically ex- 
pected behavior, the mean squared displacement (|1.1[) . 
the fractal dimension df = 2 — a/2 of the resulting trajec- 
tory, and the first passage time distribution. To simulate 
the confined motion, reflecting walls were considered at 
locations ±L for each coordinate. For instance in the ID 
case, if \x H (t)\ > L at some time t, the particle bounces 
back to the position x H (t) — 2\x H (t) — s\gx\(x H )L\. Sim- 
ilar reflecting conditions were taken into account in the 
multi-dimensional case. 



B. Fractional Langevin equation motion 

In simulating FLE motion, we follow the numerical 
method presented by Deng and Barkai [25J]. First, inte- 
grating Eq. (|2.15|) from to t, we obtain the Volterra 
integral equation for velocity field v(t) = dy(t)/dt 

v(t) = - 1 fit t'fv-'v^dt' 
(2H - l)m J a 

+v + H-x W (t), (3.1) 
m 

where v(t = 0) = vq. This stochastic integral equa- 
tion can be evaluated by the predictor-corrector algo- 
rithm presented in Ref. [57| with the FBM x H (t) inde- 
pendently obtained by the Hosking method. We calcu- 
lated y(t) = yo + Jo v (t')dt' by the trapezoidal algorithm. 
For discrete time steps, the equation of motion is given 
by 

Vn+i = 2/o + («o + v n+ i) +dh^Tvi, 

i=i 

dh 

= —(v n + v n+ i) + y n , (3.2) 

where dh is the time increment. When evaluating 
Eq. (|3.2[) . a reflecting boundary condition was considered 
in the sense that y n ->• y n - 2\y n - sign(y n )L\ if \y n \ > L. 

To show the reliability of our simulation, we compare 
the simulation result with the well-known solution for 
free space motion. In Fig. 1, we simulate the subdiffu- 
sion case with the parameters values, H = 5/8, m = 1, 
vq = 1, 2/o = 0, 7 = 10, k B T = 1, and dh = 0.01. From 
200 simulated trajectories, we obtain the ensemble av- 
eraged and time averaged mean squared displacements, 
(y 2 (t)) and 5 2 (A,T), and compare them with the exact 
solution Eq. (|2.19j) . Note that (y 2 (t)) should be identical 
to (5 2 (A, T), Eq. ([572]) . with t regarded as the lag time A 
due to the ergodicity of the FLE motion in free space [25[ . 
The deviation from the exact solution is markedly re- 
duced with decreasing time increment dh. With our cho- 
sen value dh = 0.01 the mean squared displacements ob- 
tained from simulation appear to be in good agreement 
with the theory. 
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FIG. 1. (Color online) The mean squared displacements 
(MSD) for FLE motion in free space. The ensemble av- 
eraged and time averaged MSDs, (y 2 (t)} and 8 2 (A,T), ob- 
tained from simulation are compared with the exact solution 
2i 2 £ 5 /4,3[-10r(5/4-l)f 5/4 ] given by Eq. (pT9|) . The ensem- 
ble averaged value was obtained from 200 simulated trajecto- 
ries. In the simulation, we chose the Hurst exponent H = 5/8, 
time increment dh = 0.01, particle mass m = 1, initial veloc- 
ity vo — 1, initial position yo = 0, friction coefficient 7 = 10, 
and k B T = 1. 



IV. FRACTIONAL BROWNIAN MOTION IN 
CONFINED SPACE 

We now turn to the investigation of the behavior of 
FBM under confinement, analyzing the mean squared 
displacement and potential crgodicity breaking. We then 
define the displacement correlation function, and finally 
study the influence of dimensionality 



A. Mean squared displacement 

For FBM in free space (x H (t) 2 ) can be estimated by 
the time averaged mean squared displacement Eq. (|1.2p 
via the exact relation [25l 



S 2 (A,T) )=2K H A 



2H 



(4.1) 



Here (•) denotes the ensemble average. In contrast to 
CTRW subdiffusion, in FBM this quantity is ergodic. 
However, as mentioned above, the approach to ergod- 
icity is algebraically slow, and we want to explore here 
whether boundary conditions have an impact on the er- 
godic behavior. Let us now analyze the behavior in a box 
of size 2L. 

In Fig. [2] we show typical curves for the time averaged 
mean squared displacement for Hurst exponent H = 1/4 



FIG. 2. (Color online) Time-averaged mean squared displace- 
ment (MSD) versus the lag time A for given values of L. The 
drawn line has slope 1/2, corresponding to the expected short 
lag time behavior for the used value H = 1/4 of the Hurst 
exponent. For L = 1, 2, and 3, five different trajectories each 
are drawn to be able to see whether the trajectories scatter. 
The simulation time is T = 2 17 « 1.3 • 10 5 . 



and three different interval lengths L. Regardless of the 
size of L, the confined environment does not affect the 
power law with exponent 2H for short lag times. More- 
over at long lag times we observe saturation of the curves 
to a value that depends on L. This behavior is distinct 



from that of the CTRW case where {6 2 {A,T)j shows a 

power law with slope 1 — a [2?], HH . 

One can estimate the saturated value as a function of 
L. For long A and measurement time T, the probability 
p(x) to find the particle located at x is independent of x 
due to the equilibration between the reflecting walls, and 
thus J^ L x 2 p(x)dx = L 2 /'S. The dotted lines in Fig. [2] 
represent these values. 

We observe that the scatter between different single 
trajectories becomes more pronounced when the interval 
length is increased. In fact the scatter is negligible for 
L = 1 while it is quite appreciable for L = 3, even though 
the slope of all curves at finite L converges to a horizon- 
tal slope, with an amplitude close to the predicted value 
L 2 /3. We also note that the scatter depends on the total 
measurement time T. For given L it tends to be reduced 
as we increase T. This effect will be discussed quantita- 
tively in detail using the ergodicity breaking parameter. 



B. Ergodicity breaking parameter 

In contrast to CTRW subdiffusion, FBM in free space 
is known to be ergodic (25j . The time averaged mean 
squared displacement traces displayed in Fig. [2] exhibit 
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no extreme scatter as known from the CTRW case. This 
implies that ergodicity is indeed preserved for confined 
FBM. We quantify this statement more precisely in terms 
of the ergodicity breaking parameter |23f 

Eb(A,T) = -S ; 1 , 9 , (4.2) 

S 2 (A,T) 

where limT-j-oo Eb{T) = is expected for crgodic sys- 
tems. For the case of free FBM, Deng and Barkai ana- 
lytically derived that Eb decays to zero as 



E B (A,T) 



A 

T 



for < H < 



for H = 



(4.3) 



4-4H 



for f 



<H < 1, 



for long measurement time T [25j. 

We numerically investigate the boundary effects on the 
ergodicity breaking parameter. First, in Fig. [3] we eval- 
uate Eb as function of the lag time A from 200 FBM 
simulations for each given L. The dotted line represents 
the expected free space behavior Eb ~ A, which is nicely 
fulfilled by the data at shorter times and sufficiently large 
L. At longer times or small L the results show that 
Eb behaves very differently when confinement effects arc 
present. The plateau in Eb is related to the saturation 
of the curves for the mean squared displacement, Fig. [5] 
As the motion is restricted by the walls roughly above a 
crossover lag time A cr = (L 2 /2ifff) 1 / 2ff , the ergodicity 
breaking parameter Eb levels off at A > A cr . The sharp 
increase at the end of the curve is due to the singularity 
when the lag time reaches the size of the overall mea- 
surement time T, which would disappear in the infinite 
measurement time. 

In Fig. 2] we show Eb for given A as function of 
the measurement time T for the same choice of interval 
lengths, L = 1,3,5, and 10. For short lag times A, all 
Eb curves coincide and decay as T~ x , in complete anal- 
ogy to the free space motion (dotted line). In the case of 
longer A the general trend is that Eb decays like T~ x , 
unaltered with respect to the free case. However, there is 
a sudden decrease in Eb for the smallest interval size, for 
L = 1. One can understand this behavior by observing 
the Eb curve for L = 1 in Fig. [3J as the fluctuations of 
the mean squared displacement are strongly suppressed 
due to the tight confinement in this case, Eb has almost 
no dependence on A for A cr < A < T and the saturated 
value is quite small compared to those for other cases. 
Therefore, the curves for L = 1 appear disconnected from 
the other curves. Corresponding to the approximate in- 
dependence of the L = 1 curve for A > 10 in Fig. [31 
we observe in Fig. [4] that at longer times T the L = 1 
curves approach each other. Only at T w A these curves 
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FIG. 3. (Color online) Ergodicity breaking parameter Eb ver- 
sus lag time A for L = 1, 3, 5, and 10 (from bottom to top) 
with Hurst exponent H = 1/4. The overall measurement time 
is T — 2 14 ~ 1.6 • 10 4 . The dotted line with slope 1 represents 
the theoretical expectation Eb — A in free space. For each L 
the curve was obtained from 200 single trajectories. The ver- 
tical lines show the crossover lag time A cr = (L 2 /2Kh) 1 ^ 2H 
for L = 1, 3, and 5. 



separate, as then 6?(T,T) = [xf(t + T) - xf (t)] 2 , and 
Eb is evaluated with the same small number of squared 
displacement data. Note that the splitting of the Eb 
curve can be also observed for larger L at As larger than 
A cr under longer total measurement time T as other Eb 
curves also have corresponding constant saturation val- 
ues for A > A cr (= (L 2 /2KhY/ 2H ) which increases with 
the size L. 



C. Displacement correlation function 

As explained for the stochastic properties of FBM in 
Sec. HH the position autocorrelation (x H {t\)x H (t 2 )) ex- 
plicitly depends on t\ and t 2 as well as their difference, 
|£i — *a I - It is therefore not an efficient quantity to es- 
timate directly from experimental or simulations data. 
However, the correlation function of the displacements 



Sx H (t, A) = x H (t + A) - x H (t) 



(4.4) 



depends only on the time interval A of the displacement, 

(6x H (t, A)Sx H {t - A, A)) = K H {2 2H - 2)A 2H (4.5) 

for free FBM. This relation is derived in App. [3] This 
quantity is anticorrelated for < H < 1/2 (subdiffu- 
sive motion), uncorrelated for H = 1/2 (normal Brown- 
ian motion), and positively correlated for 1/2 < H < 1 
(supcrdiffusion) . As Eq. (|4.5[) does not depend on the 
measurement time T the value of the ensemble averaged 
value is identical to the corresponding time average. 
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FIG. 4. (Color online) Ergodicity breaking parameter Eb 
versus overall measurement time T at given lag time A = 1, 
100, and 1000. The dotted line depicts a power law with slope 
-1, representing the analytic behavior Eb — T" 1 in free space. 
For each A, the Eb curves are drawn for different interval 
lengths L = 1, 3, 5, and 10. Each curve was obtained from 
200 single trajectories, and the Hurst exponent was H — 1/4. 



10 1 




10" 3 -l — — — — — -I 

1 10 100 1000 10000 

lag time A 

FIG. 5. (Color online) Absolute value of the displacement 
correlation (DC) versus lag time A for L = 1, 2, and 3 (from 
bottom to top) with a slope proportional to A 2H (dotted line) . 
Total measurement time is T = 2 17 ~ 1.3 ■ 10 5 , and the Hurst 
exponent is H — 1/4. Each curve was obtained via time 
averaging from a single particle trajectory. 

Wc present the time averaged displacement correlation 
for confined subdiffusive motion (H = 1/4) in Fig. [5] Be- 
cause of the negativity of expression (|4.5p the absolute 
value of the displacement correlation is drawn in the log- 
log representation. At short lag times the slope of the 
correlation functions is proportional to 2H as expected 



FIG. 6. (Color online) Time-averaged mean squared displace- 
ment (MSD) curves versus lag time A for L = 1 in ID, 2D, 
and 3D space (from bottom to top), with a slope A 2H (dot- 
ted line). For each dimension, 5 trajectories were drawn with 
total measurement time T = 2 17 and H = 1/4. 

from Eq. (|4.5p . However, at long lag times, we inter- 
estingly observe fluctuations of the correlations around a 
constant value, reflecting the confinement of the motion. 

D. Dimensionality 

To mimic the anomalous diffusion of particles in- 
side biological cells, we also simulate two- and three- 
dimensional FBM based on Eq. (|2.7|) in the presence of 
reflecting walls. In free space, the ensemble average of 
the time averaged mean squared displacement is simply 
given by 

(P(K^)) = ¥ ^£ A ([ x H (t + A)- x H (t)] 2 )dt, 

= 2dK H A 2H , (4.6) 

i.e., it is additive as for the ensemble average. This be- 
havior is indeed observed in Fig. [5] where five different 
mean squared displacement curves are drawn for L = 1 
in ID, 2D, and 3D, respectively. Only the height of the 
curves are affected by the dimensionality. There is no 
noticeable difference in the scatter of the curves. 

We further investigate the effects of dimensionality on 
the scatter of the mean squared displacement curves, as 
possibly the strong scatter observed in experiments 0- 
flOj may also occur for FBM in higher dimensions. To see 
the effect of dimensionality on the ergodicity behavior 
we measure Eb versus lag time for one-, two-, and three- 
dimensional embedding dimension for the same values of 
L and H. Interestingly, the result shows that Eb tends to 
decrease with increasing dimensionality d, meaning that 
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lag time A 



FIG. 7. (Color online) Rescaled ergodicity breaking param- 
eter cLEb versus lag time A for interval size L — 5 for ID, 
2D, and 3D space (from top to bottom), with Hurst exponent 
H = 1/4 and measurement time T = 2 14 . Eb was evaluated 
from 200 trajectories with different initial positions. 

for FBM big scatter is not caused by higher dimensions 
in presence of reflecting walls. In fact, from Eqs. (|4.2|) 
and (|4.6|) . we can analytically derive the relation 

EbW . M^Jl, (4 . 7) 

which still holds in the case of confined motion (see ap- 
pendix B for the derivation) . This relation is numerically 
confirmed in Fig. [7] where three Eb curves collapse upon 
rescaling by dEs(d). According to this relation, we ex- 
pect that ergodic behavior obtained in one-dimensional 
confined motion (Figs. [3] and [4j> will also be present in 
multiple dimensions with a factor of l/d. 

V. FRACTIONAL LANGEVIN EQUATION 
MOTION IN CONFINED SPACE 

In this section we analyze FLE motion under confine- 
ment. Due to the different physical basis compared to 
FBM, in particular, the occurrence of inertia, we observe 
interesting variations on the properties studied in the pre- 
vious section. 



A. Mean squared displacement 

Using the correlation function [55| 

(y(h)y(h)) = ^\tlE 2 H,3{-it?) + tlE 2 x^it?) 

-{tz-hfE^^-^-hf 1 )], (5.1) 




0.01 0.1 1 ^ 10 100 



lag time A 



FIG. 8. (Color online) Time-averaged mean squared displace- 
ment (MSD) versus lag time A. The two dashed lines repre- 
sent the two asymptotic scaling behaviors ^<5 2 (A,T)^ ~ A 2 

and (d 2 (A, T)\ ~ A 2 ' 2 " . For each given L = 1/2, 1, and 3, 

five different trajectories are drawn to visualize the scatter. 
As a reference curve for motion in free space, the curve for 
L — 100 (line) is drawn. In the simulation, we chose the Hurst 
exponent H = 5/8 [NB: for the FLE this means subdiffusion] , 
time increment dh = 0.01, particle mass m = 1, initial veloc- 
ity vo = 1, initial position yo — 0, friction coefficient 7 = 10, 
and k B T = 1. 

one can show analytically that, similarly to the FBM, the 
ensemble averaged second moment (y 2 {t)) is identical to 

its time averaged analog ^<S 2 (A)^}, for all A in free space, 

namely 

(gSfcTf) = 2^A 2 ^ 3 (- 7 A^). (5.2) 

Thus, the time averaged mean squared displacement 
turns over from a ballistic motion 

^ 2 (A,T)^ ~ A 2 (5.3) 

at short lag time to the subdiffusivc behavior 

^ 2 (A,T)J) ~ A 2 - 271 (5.4) 

at long lag times, in free space. 

We numerically study how this scaling behavior is af- 
fected by the confinement. Figure [8] shows typical curves 
for the time averaged mean squared displacement, for in- 
terval sizes L = 1/2, 1, 3, and 100 (regarded as free space 
motion) with identical initial conditions and Hurst expo- 
nent H = 5/8. The results are summarized as follows: 

(1) We observe both scaling behaviors, ( S 2 (A,T)) ~ A 2 
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turning over to ~ /\2-2H ^ f or connnec i FLE motions. 
(2) For narrow intervals, the curves eventually reach the 
saturation plateau within the chosen total measurement 
time T. The saturation values are approximately L 2 /3. 
For interval size L = 1/2, the saturated value is notice- 
ably larger than L 2 /3, which appears to be caused by 
multiple reflection events on the walls. The same behav- 
ior is observed in the FBM case when considering a large 
value of H > 1/2, or very narrow intervals for the given 
H = 1/4. (3) As in the case of FBM, the scatter becomes 
pronounced as the interval length increases. 



B. Ergodicity breaking parameter 

From a simple argument and simulations it was shown 
in Ref. [2f| that the FLE and FBM mean squared 

displacements are asymptotically equal, (S 2 (y) 



S 2 (x H )), similarly for the ergodicity breaking param- 
eter, E B {y) ~ E B (x H ). Here the asymptotic equivalence 
is valid at long measurement times T, and the derivation 
holds for motion in free space. From 200 trajectories of 
the mean squared displacement we measure the ergodic- 
ity breaking parameter Eb as function of lag time A for 
interval lengths L = 1/2, 1, 3, and 100 in Fig. M The 
behavior is similar to the corresponding curves for FBM, 
displayed in Fig. [3j Eb significantly deviates from the 
reference curve for free space motion (i.e., the longer A 
behavior for L = 100 and the drawn power law ~ A), 
due to the confinement effect. Eb tends to decrease with 
smaller L for the same value of A. However, the plateau 
at short lag times that is still observed for L = 100 (re- 
garded as free space motion) , is due to the initial ballistic 
motion of FLE. In that regime the initially directed mo- 
tion renders the random noise effect negligible. 

C. Displacement correlation function 



Using the correlation function {y(ti)y(t2)}, we analyti- 
cally obtain the displacement correlation function in free 
space in the form (refer to App.[X]for the derivation) 

.k R T 



(5y{t, A)5y(t - A, A)) = 4^-A^ 27? 3 [- 7 (2A)^] 

A 2 ^ r _ a 2~H 



so that we observe the following asymptotic behavior 
(5y{t,A)5y(t-A,A)) 



(5.5) 



2k B T A2 



mr(3) 



for A -> 0, 



ra7r(3 - 2H) 



(5.6) 



where 8y(t,A) = y(t + A) — y(t). Above expression 
shows that the displacement correlation has two distinct 
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FIG. 9. (Color online) Ergodicity breaking parameter Eb 
versus lag time A for interval sizes L = 1/2, 1, 3, and 100. The 
straight line corresponds to the free space behavior Eb — A. 
Each curve was obtained from 200 single trajectories, with 
the same parameter values used in Fig. [8] 



scaling behaviors. At short lag times, it grows like A 2 
and is positive, due to the ballistic motion. At long lag 
times, it is negative in the domain 1/2 < H < 1, ex- 
hibiting the same subdiffusive behavior as observed for 
FBM [cf. Eq. ([Q]) ] when we replace H -t 2 - 2H. Note 
that to bridge these two scaling behaviors the displace- 
ment correlation passes the zero axis at A = A c that 



satisfies IE. 



2H,3 



-7(2A C ) 



2 if 



E. 



2H,3 



-7A 2 



2 if 



space. For small 7, we find approximately 



A, 



T(2H + 3) 



2T(2H - 1)(2 2H + 1 - 1) 7 



1/2H 



in free 



(5.7) 



such that it becomes exactly the momentum relaxation 
time 771/7 f° r normal Brownian motion (H = 1 /2). In the 
limit H — > 1, A c goes to infinity to satisfy the equality 



2E. 



2 if, 3 



-7(2A ( 



\2ff 



= E. 



2 if, 3 



- 7 A 



2 n 



Thus, A r can 



be interpreted as the typical timescale for the persistence 
of the ballistic motion. 

Figure [TU] shows (a) the displacement correlation ver- 
sus lag time A, and (b) the absolute value of the dis- 
placement correlation as function of A, for L = 1/2, 1, 
3, and 100. The scaling properties derived in Eq. (|5.6j) 
are indeed observed. At short lag times all curves are 
positive and scale like ~ A 2 , before decreasing to zero. 
In the long lag time regime the displacement correla- 
tion becomes negative and the predicted scaling behavior 
~ A 2 ~ 2H is observed. For small intervals (L — 1/2 and 
1), it is saturated due to the confinement effect as seen 
in the case of FBM. 
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FIG. 10. (Color online) (a) Displacement correlation (DC) 
versus lag time A for L — 1/2, 1, 3, and 100 (from bottom 
to top), (b) Absolute value of the displacement correlation 
as a function of A for the given values of L. The two slopes 
correspond to the limiting behavior A 2 and A 2_2fl \ In (a) 
and (b), the ensemble averaged curves were obtained from 
200 different time averaged displacement correlation curves. 
Same parameter values as in Fig. [8] 



D. Dimensionality 



FBM and FLE motions in confined space. In particular 
we analyzed the effects of confinement and dimensional- 
ity on the stochastic and ergodic properties of the two 
processes. Interestingly for both stochastic models, the 
confinement tends to decrease the value of the ergodicity 
breaking parameter Eb compared to that in free space. 
The same trend is observed for increasing dimensionality. 
Correspondingly the scatter of time averaged quantities 
between individual trajectories is quite small, apart from 
regimes when the lag time A becomes close to the over- 
all measurement time T and the sampling statistics for 
the corresponding time average become poor. The relax- 
ation of the ergodicity breaking parameter as function of 
measurement time is quite similar to previous results in 
free space. Wc conclude that neither confinement nor di- 
mensionality effects lead to the appearance of significant 
ergodicity breaking or scatter between single trajectories. 

The displacement correlation function introduced here 
is a useful quantity that can be easily obtained from sin- 
gle particle trajectories. It can be used as a tool to dis- 
criminate one stochastic model from another. For subd- 
iffusive motion governed by FBM and FLE motion, the 
displacement correlation should be negative and saturate 
in the long measurement time limit due to the confine- 
ment. Notably, the negative decrease (~ — A Q ) with lag 
time A and anomalous diffusion exponent a is an intrin- 
sic property of FBM and FLE displacement correlations 
which is clearly distinguished from that of CTRW subdif- 
fusion. In the latter case, the subdiffusive motion occurs 
due to the long waiting time distribution between suc- 
cessive jumps and there is no spatial correlation between 
them, so that displacement correlation only fluctuates 
around zero with time. FLE motion can be distinguished 
from FBM motion since the displacement correlation has 
a positive value at short times due to the ballistic motion 
in the FLE model. 



In the case when the memory tensor is diagonalized, 
each coordinate motion is independent and FLE mo- 
tion exhibits qualitatively the same behavior as shown 
in the case of FBM with increasing dimensionality From 
the mean squared displacement curves, the same scal- 
ing behavior is expected with more elevated amplitude 
for higher dimensionality. In fact, when each coordinate 
motion is decoupled, a d-dimensional motion effectively 
increases the number of single trajectories d times com- 
pared to the one-dimensional case. Therefore, the scatter 
in the mean square displacement curves decreases with 
increasing dimensionality and the ergodicity breaking pa- 
rameter Eb is expected to follow the relation Eq. (|4. 7[) . 

VI. CONCLUSION 

Motivated by recent single particle tracking experi- 
ments in biological cells, in which confinement due to 
the rather small cell size becomes relevant, we studied 
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Appendix A: Derivation of the displacement 
correlation function 

In this appendix we derive analytical expressions for 
the displacement correlations, Eqs. (|4.5p and (|5.5|) . For 
a stochastic variable x, we define 

Sx(t, A) = x(t + A) - x(t). (Al) 

The displacement correlation is then given by 

(Sx(t, A)Sx(t — A, A)) = 

(x(t + 2A)x(t + A)) - (x(t + 2A)x(t)) 
-(x(t + A) 2 ) + (x(t + A)x(t)). (A2) 
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We now calculate this expression for FBM and FLE mo- Expanding the generalized Mittag-Lefflcr function 

E 2l j 3 (x) « l/r(3) + x/T(2H + 3) + • • • for x < 1 the 
displacement correlation is approximated as 



tions. 



FBM 



(5x(t, A)5x(t - A, A)) 



2 k B T_ A2 



T(3) m 



(A7) 



For FBM (x(f) = x H (t)), we use the expression 



at short lag times. With the expansion E 2 jj 3 (— x) 



(x(fi)x(f 2 )) = K H {t\ H +t% H - |f 2 - h\ ztl ) (A3) l/xr(3 - 2H) for x > 1 the long lag time behavior of 



\2H 



for the autocorrelation. With this we readily obtain the 
result 



the displacement correlation is obtained as 

2 2 ~ 27? -2 k B T 



(6x{t,A)6x(t-A,A)) 



(6x{t, A)Sx(t - A, A)) = K H {2 ZH - 2)A 2tl . (A4) 
2. FLE 

For FLE (x(t) = y(t)), we use the correlation function 



7 r(3-2if) m 



A 



2-2H 



(A8) 



Notejhat the prefactor (2 2 ~ 2ff - 2)/r(3 - 2H) is zero 
for H — 1/2 and then becomes increasingly negative, 
saturating at the value 1 for H = 1. 



Appendix B: Derivation of Equation (|4.7I 



(x(fx)x(f 2 )) 



~ [tlE 2 H* (-jt 2H ) +tlE 2 - n J- 1 tl H ) 



From the definition of the time averaged mean squared 



-(f 2 -hfE^^-^h -h\ 2H )]. (A5) displacement, Eq. gU), we expand ^((5 2 (A,T) 

The displacement correlation is then obtained as the form 



(6x(t, A)Sx(t - A, A)) 
= ^A*E- 



2H,3 



- 7 (2A) 



2H 



—^-A 2 E. 



2H,3t 



-j A 



2/fi 



(A6) 



<5 2 (A,T) 



i-T—A i-T—A . d, 



(T- 



l 7AjL I {X>^ + A) - (*i)] 2 K H (*2 + A) - xf (t 2 )] 2 ) 

+ £>f (*i + A ) - (*i)] 2 )(K ff (*2 + A) - xf (t 2 )] 2 )}dhdt 2 . (Bl) 



Using the Isserlis theorem for Gaussian process with zero mean [48(: 

(x(fx)x(f 2 )x(f 3 )x(f 4 )) = (a(ti)x(ta))(aj(t 3 )a:(t4)> + (x(fi)x(f 3 )>(x(f 2 )x(f 4 )) + (x(fi)x(f 4 )}(x(f 2 )x(f 3 )), 
the first term in the braces in Eq. (|B1[) can be rewritten as 



(B2) 



£>f (f x + A) - (f0] 2 [xf (f a + A) - xf (i 2 )] 2 > = £>f (fx + A) - (fr)] 2 )([xf (t 2 + A) - xf (f 2 )] 2 > 
i=i i=i 

+2^([xf (f x + A) - xf (fx)][xf (f 2 + A) - xf (i 2 )]> 2 . 



(B3) 



In this expression, we note that the sum of the second term in Eq. (|B1[) and the first term in Eq. (|B3|) yields 
<* 2 (A,T)) 2 : 



<5 2 (A,T) 



(r-A) 2 7 



T-A i>T— A 



([x H (f x + A) - x ff (fx)] 2 )([x ff (f 2 + A) - x H (t 2 )] 2 )dMf 2 



Di 



(B4) 
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where we used the property 

X>f + A) - xf( tl )] 2 )([xf(t 2 + A) - xf {t 2 )f) = d 2 {[x H {h + A) - x H ( tl )] 2 )([x H (t 2 + A) - x H (t 2 )f) (B5) 



due to the independence of the motion in each coordinate direction. We also note that the expression ( I <5 2 (A, T) 



S 2 {A 1 T)j simplifies to 

((^T) 2 ) - = (r _ 2 A)2 E/ / ([af (ti + A) -af (ti)] x [xf (t 2 + A) - xf {t 2 )]) 2 dt x dt 2 

= d[{{pf) lD mi D }. (B6) 

From Eqs. ()B4[) and (|B6J) . the ergodicity breaking parameter follows the general relation 
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